refract_time=100;
l_win=1000;
r_win=5000;
min_dur=50;


mo_trig_saccades=[];
active_cells=[];
cell_cnt=0;

rec_dur=[];

do_layer=2;

saccades_per_sec=nan(length(AV_data),5);
saccades_total_time=nan(length(AV_data),5);
saccades_number=nan(length(AV_data),4);
running_fract=nan(length(AV_data),1);

running_speed_hist=nan(84,length(AV_data));
start_running_after_pb_onset=[];
vel_snip=nan(15001,2500);
vel_snip_cnt=0;

for ind=1:length(AV_data)
    disp(num2str(ind))
    if AV_data(ind).layer==do_layer & AV_data(ind).ACode==1
        pb_onsets=[];
        if length(AV_data(ind).shutter_open)==2
            pb_onsets(1)=AV_data(ind).shutter_open(1);
            pb_onsets(2:4)=find(diff(AV_data(ind).PS(AV_data(ind).shutter_open(2):end)>0.1)>0)+AV_data(ind).shutter_open(2);
            pb_offsets=find(diff(AV_data(ind).PS(AV_data(ind).shutter_open(2)+1000:end)>0.1)<0)+AV_data(ind).shutter_open(2);
            pb_dur=pb_offsets(1)-pb_onsets(2);
            interim_runnig(3)=sum(AV_data(ind).velP_ind(round(pb_offsets(1)/10):round(pb_onsets(3)/10)));
            interim_runnig(4)=sum(AV_data(ind).velP_ind(round(pb_offsets(2)/10):round(pb_onsets(4)/10)));
        elseif length(AV_data(ind).shutter_open)==3
            pb_onsets(1)=AV_data(ind).shutter_open(1);
            pb_onsets(2:4)=find(diff(AV_data(ind).PS(AV_data(ind).shutter_open(3):end)>0.1)>0)+AV_data(ind).shutter_open(3);
            pb_offsets=find(diff(AV_data(ind).PS(AV_data(ind).shutter_open(3)+1000:end)>0.1)<0)+AV_data(ind).shutter_open(3);
            pb_dur=pb_offsets(1)-pb_onsets(2);
        elseif length(AV_data(ind).shutter_open)>=4
            pb_onsets=AV_data(ind).shutter_open+1000;
            pb_dur=pb_onsets(2)-pb_onsets(1)-2000;
        else
            pb_onsets=[];
        end
        pb_onsets=pb_onsets+1000;
        pb_dur=pb_dur-2000;
        
        rec_dur(ind)=AV_data(ind).frame_times(end)-AV_data(ind).frame_times(1);
        
        type=(AV_data(ind).velP_ind&AV_data(ind).velM_ind)+2*(AV_data(ind).velP_ind&~AV_data(ind).velM_ind)+3*(~AV_data(ind).velP_ind&AV_data(ind).velM_ind);
%         type_offsets=find(diff(type>0&type<3)==-1);
%         for knd=1:length(type_offsets)
%             type(max(1,type_offsets(knd)-100):type_offsets(knd)+100)=0;
%         end
        
        for knd=1:size(AV_data(ind).raw_act,2)
            if max(AV_data(ind).raw_act(:,knd))<100
                cell_cnt=cell_cnt+1;
                std_noise=sqrt(sum((AV_data(ind).raw_act(AV_data(ind).raw_act(:,knd)<1,knd)-1).^2/sum(AV_data(ind).raw_act(:,knd)<1)));
                std_noise_all(cell_cnt)=std_noise;
                
                for jnd=1:4
                    act_diff_tmp=diff([0 (AV_data(ind).raw_act(type(round(AV_data(ind).frame_times/10))==jnd-1,knd)<(1-3.72*std_noise))' 0]);
                    if ~isempty(max(find(act_diff_tmp==-1)-find(act_diff_tmp==1))>0)
                        active_cells(cell_cnt,jnd)=max(find(act_diff_tmp==-1)-find(act_diff_tmp==1))>8;
                    else
                        active_cells(cell_cnt,jnd)=0;
                    end
                end
                %                 active_cells(cell_cnt)=sum(AV_data(ind).raw_act(:,knd)<(1-4.26*std_noise))>20;
            end
        end
        
        [mo_start]=parse_cond_starts(AV_data(ind).velP_ind,min_dur,refract_time);
        %[vs_start]=parse_cond_starts(~AV_data(ind).velP_ind&AV_data(ind).velM_ind,min_dur,refract_time);
        [ps_raw]=find_ps_during_running(AV_data(ind).velP_ind,AV_data(ind).velM_ind);
        [ps_start]=parse_cond_starts(ps_raw,min_dur,refract_time);
        
        jjj(ind)=length(ps_start);
        
        % this will look for pb offsets during non-running
        [vs_raw]=find_pb_during_sitting(AV_data(ind).velP_ind,AV_data(ind).velM_ind);
        [vs_start]=parse_cond_starts(vs_raw,min_dur,refract_time);
        
        for jnd=1:length(vs_start)
            try
                if sum(AV_data(ind).velP_ind(round(vs_start(jnd)/10)-500:round(vs_start(jnd)/10)))<300
                    vel_snip_cnt=vel_snip_cnt+1;
                    vel_snip(:,vel_snip_cnt)=AV_data(ind).velP(vs_start(jnd)-5000:vs_start(jnd)+10000);
                    starts_within_range=abs(vs_start(jnd)-mo_start)<20000;
                    start_running_after_pb_onset(end+1:end+sum(starts_within_range))=mo_start(starts_within_range)-vs_start(jnd);
                end
            end
        end
        
        running_fract(ind)=sum(type>0&type<3)/length(type);
        
        pb_onsets_m=round([pb_onsets length(AV_data(ind).PS)]/10);
        for knd=1:length(pb_onsets_m)-1
            type_session=type(pb_onsets_m(knd):pb_onsets_m(knd+1));
            running_fract_session(knd,ind)=sum(type_session>0&type_session<3)/length(type_session);
        end
        
        running_speed_hist(:,ind)=hist(-AV_data(ind).velP,[0:0.0012:0.1]);xlim([0 0.1]);
        running_speed_hist(:,ind)=running_speed_hist(:,ind)/length(AV_data(ind).velP);
        
        
        saccade_times_tmp=AV_data(ind).iframe_times(AV_data(ind).saccade_frames);
        
        for jnd=0:3
            saccades_per_sec(ind,jnd+1)=sum(type(round(saccade_times_tmp/10))==jnd)/sum(type==jnd)*100;
            saccades_total_time(ind,jnd+1)=sum(type==jnd)/100;
            saccades_number(ind,jnd+1)=sum(type(round(saccade_times_tmp/10))==jnd);
        end
        saccades_per_sec(ind,5)=length(saccade_times_tmp)/AV_data(ind).iframe_times(end)*1000;
        saccades_total_time(ind,5)=AV_data(ind).iframe_times(end)/1000;
        
        use_saccades=ones(length(saccade_times_tmp),1);
        for jnd=1:length(saccade_times_tmp)
            if min(abs(saccade_times_tmp(jnd)-pb_onsets+1000))<1000
                use_saccades(jnd)=0;
            end
        end
        saccade_times_tmp=saccade_times_tmp(logical(use_saccades));
        
        for jnd=1:length(ps_start)
            use_saccades=find(abs(ps_start(jnd)-saccade_times_tmp)<5000);
            mo_trig_saccades(end+1:end+length(use_saccades))=saccade_times_tmp(use_saccades)-ps_start(jnd);
        end
    end
end

% saccade detection problems
saccades_per_sec(3,:)=NaN;
saccades_per_sec(10,:)=NaN;
saccades_total_time(3,:)=NaN;
saccades_total_time(10,:)=NaN;
saccades_number(3,:)=NaN;
saccades_number(10,:)=NaN;

nanmean(saccades_per_sec)
nanstd(saccades_per_sec)/sum(~isnan(saccades_per_sec(:,5)))

figure;
bar([0:0.0012:0.1],nanmean(running_speed_hist'),1)
set(gca,'xtick',[0:0.012:0.1],'xticklabel',[0:10:83])
xlabel('Running speed [cm/s]')
ylabel('Fraction of time')
xlim([-0.001 0.1])
ylim([0 0.035])
ylim([0 0.9])



figure;
hist(mo_trig_saccades,100)
box off
title('PS triggered saccades')
ylabel('Number of saccades')
